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We study lattice QCD with two flavors of dynamical domain wall quarks. With renormalization group motivated 
actions, we find chiral symmetry can be preserved to a high degree at lattice cut off of a -1 ~ 2 GeV even for fifth 
dimension size as small as L s — 12. In addition two new steps are introduced to improve the performance of the 
hybrid Monte Carlo simulation. 



1. Introduction 

An obvious source of systematic error in nu- 
merical lattice QCD at present is the quenched 
approximation. Use of domain wall fermion 
(DWF) [Q dynamical quarks to lift this approxi- 
mation is attractive because DWF preserve the 
flavor and chiral symmetries both on and off 
shell on the lattice before the continuum limit 
is taken. Evidently successful implementation 
would bring many advantages, such as detailed 
study of flavor-singlet meson spectrum and its 
relation to the U(1)a anomaly. However some 
recent exploratory calculations Q on coarse lat- 
tices with cut off of a^ 1 < 1 GeV show the chiral 
symmetry is rather poorly realized: impractically 
large L s ~ 100 are necessary to bring the resid- 
ual symmetry violation under control. In contrast 
in this work we go to a higher cutoff of a -1 ~ 2 
GeV with renormalization-group (RG) motivated 
gauge actions and find that the residual violation 
can be made sufficiently small even at L s = 12. 

To cope with considerably larger demand for 
computing resources brought by this new ap- 
proach, we introduce two improvements in the 
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hybrid Monte Carlo (HMC) algorithm which ac- 
celerate the simulation reported below on the 
QCDSP by about a factor of three. 

2. New algorithms 

The lattice QCD partition function with two 
flavors (Np =2) of DWF quarks is 



vu det(D^D)\ m=mf c _ Sa 



(1) 



det{DW)\ n 

where D is the 5d Dirac operator with 4d bare 
fermion mass m and S g the gauge action. The 
determinant in the denominator in Eq. [l] comes 
from the Bosonic Pauli-Villars regulators, <f)pv, 
whose mass is set to 1 to cancel the dominating 
contribution from the fermions on the bulk 5d 
lattice. This allows us to study 4d QCD with 2 
flavors of light quarks localized at the boundary 
and coupled to the 4d gauge field. The original 
fermion action employed in pJ2j was 



Sf = Y^i-M^D 



(2) 



t> PV (tfD)\ 



m=l(PPV 
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where <$>f is the pseudo- fermion field. By inte- 
grating out 4>f and 4>pv, we obtain the partition 
function Eq. |l|. We chose hybrid Monte Carlo phi 
algorithm (HMC-</>) Q to generate the Markov 
chain of gauge configurations. 

Noticing equivalence between ^(bIb) an< ^ 
{det[B(A Jf A)- 1 B^}}- 1 , we could also obtain the 
partition function Eq. [l] from a new action ^ 

S' = -4> F [D\ m =i{D ^D)\ m l mf D\l^]<f>F (3) 

4 This part of the work is done mainly by C. Dawson. 
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where we need only one bosonic field <f>F- 
Even/odd preconditioning is implemented in ei- 
ther case. 

The benefits of the new action and the asso- 
ciated force term are first the unwanted contri- 
bution from 5d bulk fermion is now canceled ex- 
plicitly while it was canceled only stochastically 
in the original force term. This leads to better 
acceptance in the Metropolis accept/reject step, 
which increased by about 20% in a realistic sim- 
ulation. Secondly, the source vector, <p F , of the 
Dirac equation solved in the force term calcula- 
tion is modified to D^\ m —i(pp as seen in Eq. ^, 
and it turns out that the number of conjugate 
gradient (CG) iterations is 20 - 30 % smaller than 
that of the original action. 

The second improvement || is to forecast the 
solution of the CG at each step of a HMD trajec- 
tory from the solutions of previous steps. In each 
of trajectories, the system obeys the smooth clas- 
sical molecular dynamics (MD), so future solu- 
tions can be estimated from past ones. We choose 
the Minimal Residual Extrapolation method of 
H with a repeated Gram-Schmidt orthogonaliza- 
tion to ensure orthogonality. The solution is fore- 
casted by minimizing the norm of the residual 
vector within the subspace spanned by previous 
solutions. Then the forecasted solution is passed 
into the ordinary CG as an initial input vector. 

The manner of convergence of the CG based 
on different numbers of previous solutions in the 
forecasting is shown in Fig. The squared norm 
of the residual vector of a solution (vertical axis) 
is already small, 0(1), at zero CG count (hori- 
zontal axis) when the number of previous solution 
used in the forecast (N p ) is larger than one. The 
CG iterations required for convergence decreases 
monotonically with increasing N p , up to 7. 

The reversibility required to prove the detailed 
balance in the Markov chain is checked by flip- 
ping the sign of the conjugate momentum field 
at the end of a trajectory and evolving the con- 
figuration backward by the same number of MD 
steps. The resulting configuration obtained at the 
end of the reversed process differs from the initial 
configuration by less than 2 x 10~ 5 in elements of 
SU(3) link variables for a CG convergence crite- 
ria, |res|/|src| < 10~ 8 . By this trick with N p = 7 
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Figure 1. The squared norm of the residual vec- 
tor as a function of CG steps. The initial solution 
vector is precognized based on previous N p solu- 
tions. 

the number of total CG count is 55 % of the orig- 
inal code (using N p — 1) and about 35 % of a 
simulation without the forecasting. 

3. Physical scale of the dynamical lattice 

Observing the large chiral symmetry violation 
for coarser lattice spacing, a -1 < 1 GeV j2|, we 
decided to perform a large scale simulation with 
target a" 1 ~ 2 GeV. To quickly determine pa- 
rameters for this target lattice spacing, we ex- 
amined the parameter space on an 8 3 x 4 lat- 
tice and compared the meson spectrum to that 
obtained with the dynamical staggered fermions 
on the same lattice size assuming that scaling vi- 
olations from the two fermion actions are mild. 
We also explored the dependence of the residual 
chiral symmetry breaking with various RG im- 
proved gauge actions. We found at fixed a -1 ~ 2 
GeV that chiral symmetry was better realized 
for larger negative values of the coefficient of the 
(1 x 2) rectangular plaquette although the depen- 
dence of the residual mass, m res , was much more 
mild than for quenched simulations 

For our most intensive simulation, we chose 
DBW2 gauge action at = 0.80 with N F = 2 
on a 16 3 x 32 lattice. The sea quarks have mass 
am^ vn ^ = 0.02, domain wall height M5 = 1.80, 
and L s = 12. An trajectory consists of 50 leap- 
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frog steps and advances the MD time, r, by 0.5. 

The average number of CG iterations to fulfill 
the convergence criteria of |res|/|src| < 10 -8 is 
550 for N p = 1 and 314 for N p = 7. The accep- 
tance ratio of the accept/reject step is about 77(3) 
%. On a 4096 node partition (~ 205 GFLOPS 
peak) of the QCDSP, it takes about 1 hour per 
trajectory. 

To estimate the lattice spacing within the 
statistics accumulated so far, we measure the 
heavy quark potential and the meson spectrum 
from the last 50 lattices which are separated by 
five trajectories after dropping the first 500 tra- 
jectories. We monitor the chiral condensate (qq) 
as a function of the trajectory number and pick 
lattices for measurement after (qq) is considered 
to have reached its thermalized value. However, 
we note that it is not clear whether the evolution 
is thermalized at this point. 

Measuring the correlation between path or- 
dered products of gauge-fixed links in the time 
direction of length T at a space point x, Lt(x), 
we obtain the heavy quark potential e^ v ^ T = 
(L 1 T (x)L T (x + r)) 0. As shown in Fig. §, the fit 
to the data for 4 < T < 6 and y/2 < r < 7 
to the function V(r) — C — a/r + err gives 
= 0.28(2), or a" 1 w 1.7(2) GeV, and the 
Sommer scale, ro/a = 4.1(3). 

Meson propagators are calculated using valence 
quarks with L s = 12, M 5 = 1.8 and am^ a ^ = 
0.015,0.020,0.025,0.030. Extrapolating the vec- 
tor meson mass linearly, we find aMyfmJ" 1 '' — > 
0) =0.34(9) which is consistent with the results 
of the heavy quark potential within (large) er- 
ror. The squared mass of the pseudo-scalar me- 
son is well fit to the PCAC relation. Taking the 
Kaon mass from experiment and the lattice spac- 
ing from the heavy quark potential yields a sea 
quark mass approximately half the strange quark 
mass. 

Finally, the residual chiral symmetry breaking 
as determined from the pseudo-scalar mid-point 
correlator JtJ turns out to be am res (rrif — > 0) = 
1.47(7) x 10~ 3 , which corresponds to a few MeV. 



5 This calculation is done using the MILC code, version 6 
for which we thank MILC collaboration. 
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Figure 2. The heavy quark potential obtained 
from T = 4 — 5 (solid error bars) and T = 5 — 6 
(dotted error bars) as a function of distance be- 
tween quark and anti-quark with its fitted curve. 

4. Conclusion and Discussion 

In summary, we have begun investigating lat- 
tice QCD with two flavors of dynamical domain 
wall quarks. Assuming that the lattices have 
reached thermal equilibrium, we could interpret 
the configurations as QCD at roughly 0.1 fm lat- 
tice spacing, sea quarks with mass about m s /2, 
and residual chiral symmetry breaking of only a 
few MeV. This level of symmetry breaking is now 
commensurate to that observed in quenched stud- 
ies. 

Two new algorithms were introduced to en- 
hance the performance of the simulation. The 
speed up for the parameters discussed in this pi- 
lot study is about a factor of three. With this 
acceleration, further developments of theoretical 
and numerical techniques, and more computing 
horse power ||, we may reach our goal to study 
physics in the continuum limit of QCD. 
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